### =========================================================================
### Preparations
### =========================================================================

### Indicate the working path
wkpath <- ""
setwd(wkpath)


### Please indicate the path to the environmental layers
SDM_env_path <- ""
### The 10 environmental layers can be downloaded from:
### https://www.envidat.ch/dataset/sdm-env-layers-gdplants


### Indicate the path to the input occurrence files (in .csv format)
### The occurrence file should contain two columns as "x" and "y" as longitude and latitude
input_dir <- "example/cleaning_cc_occurrences"

### Indicate the path for the output maps
output_dir <- "3.2_mapping_SDM"
if (!dir.exists(output_dir)) {dir.create(output_dir, recursive = T)}

cmon.files <- list.files("./functions/sdm_functions/",full.names = T, recursive = T)
sapply(cmon.files,source)

### Please choose the strategy to generate pseudo-absence points. 
### Strategies can be chosen from:
### 'target.group','geographic', 'density', 'random',  'geo.strat', 'env.strat', 'env.semi.strat'
pa.strategy <- "geo.strat"

### If you use "target.group" strategy, please set the folder sample pseudo-absence points sampling. 
### Default is the same as the input folder.
target.group_dir <- list.dirs("2.2_cleaning_cc/cleaning_cc_occurrences", recursive = F)


### =========================================================================
### Load and install packages
### =========================================================================

#install packages globally

# set a local mirror server
options(repos=structure(c(CRAN="https://stat.ethz.ch/CRAN/")))

# Packages to load
packages <- c("raster","rgdal","maptools", "tools","dismo","cluster","class",
              "gam","gbm","randomForest","ROCR","parallel","spatstat","rgeos",
              "devtools")

# Install missing ones and load all packages
for (p in packages) {
  if(!library(package = p, logical.return = TRUE, character.only = TRUE)){
    install.packages(p)
    library(package = p, character.only = TRUE)
  } else {   
    library(package = p, character.only = TRUE) 
  }
}

# Install the if not available GitHub
if(!library(package = "PseuAbs", logical.return = TRUE, character.only = TRUE)){
  install_github("filbe87/PseuAbs")
} else {
  library(PseuAbs)
}


# Load environmental layers

  env.stk= stack(file.path(SDM_env_path, 'CHELSA_bio10_01_preped.tif'),
                 file.path(SDM_env_path, 'aridity_tranformed.tif'),
                 file.path(SDM_env_path, 'ORCDRC_M_sl1_1km_ll_tranformed.tif'),
                 file.path(SDM_env_path, 'CHELSA_FCF_1979-2013_V1.2_preped.tif'),
                 file.path(SDM_env_path, 'CHELSA_bio10_17_tranformed.tif'),
                 file.path(SDM_env_path, 'PHIHOX_M_sl1_1km_ll_preped.tif'),
                 file.path(SDM_env_path, 'CHELSA_bio10_02_preped.tif'),
                 file.path(SDM_env_path, 'CHELSA_bio10_15_preped.tif'),
                 file.path(SDM_env_path, 'CLYPPT_M_sl1_1km_ll_preped.tif'))
  
  pred_sdm <- c('Mean_Temp', 'Aridity', 'Organic_content_soil',"Frost_Change_Frequency","Prec_Driest_Quarter",
                "Soil_pH","Mean_Diurnal_Range","Prec_Seasonality","Clay_content_soil") 
  names(env.stk) <- pred_sdm
  

### =========================================================================
### Definitions
### =========================================================================

proj <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")

### =========================================================================
### Loop over species and create maps
### =========================================================================

pathfile <- list.files(input_dir, pattern = "\\.csv", recursive = T, full.names = T)
Spname <- gsub("\\.csv", "", list.files(input_dir, pattern = "\\.csv", recursive = T))

#############


for (id in 1:length(Spname)){
  
  ### Get name
  env.stk_i=env.stk
  spi_name <- Spname[id]
  print(paste("Working with",id,spi_name))

  ### Get occurrences
  
  spp_occ <- read.csv(pathfile[id])
  if (nrow(spp_occ) >= 20){
    
    spp_occ <- spp_occ[,c("x","y")]
    spp_occ <- SpatialPoints(spp_occ, proj4string=proj)
    
    # Create pseudoabsences
    pseu.abs_i <- wsl.samplePseuAbs(type=pa.strategy,
                                    n=10000, #min 1000
                                    env.stack=env.stk,
                                    add.strat=0.2,
                                    env.strat_path = file.path(output_dir, "envpath"),
                                    pres = spp_occ,
                                    template_dir= file.path(output_dir, "template_dir"),
                                    target.group_dir = target.group_dir,
                                    taxon=spi_name,
                                    force_spat_thin = "both",
                                    geores_fact=10,
                                    limdist=NA)
    print(paste("Create pseudoabsences for", spi_name, "with", PAm, "strategy."))
    
    ### =========================================================================
    ### Decide on predictor number depending on available presence observations
    ### =========================================================================

    ssize=floor(length(which(pseu.abs_i@pa==1))/10)
    
    # If there are not at least 20 presence observations, do nothing
    if(ssize<2){
      next
      # If there are less than 90 observations reduce predictor set so that
      # at least 10 observations are available per predictor
    } else if(ssize<9){
      env.stk_i=env.stk[[1:ssize]]
      pseu.abs_i@env_vars=pseu.abs_i@env_vars[, 1:ssize,drop=F]
    }
    
    ### =========================================================================
    ### Prepare data for modelling
    ### =========================================================================
    
    wip=which(pseu.abs_i@pa==1)
    blkp=sample(rep(1:3,each=ceiling(length(wip)/3)),ceiling(length(wip)/3)*3)[1:length(wip)]
    wia=which(pseu.abs_i@pa==0)
    blka=sample(1:3,size = length(wia),replace = T)
    
    blks=rep(NA,length(pseu.abs_i@pa))
    blks[wip]=blkp
    blks[wia]=blka
    
    ### =========================================================================
    ### Define SDMs
    ### =========================================================================
    
    ### Define formulas
    # GLMs
    form.glm.s=as.formula(paste("Presence~",paste(names(env.stk_i),collapse="+"))) # simple
    form.glm.i=as.formula(paste("Presence~",paste(paste0("poly(",names(env.stk_i),",2)"),collapse="+"))) # intermediate
    
    # GLM complex
    if(nlayers(env.stk_i)>1){
      cmbs<-combn(names(env.stk_i),2)
      pst<-apply(cmbs,2,paste,collapse=":")
      int.part=paste(pst,collapse="+")
      form.glm.c=as.formula(paste(paste("Presence~",paste(paste0("poly(",names(env.stk_i),",3)"),collapse="+")),int.part,sep="+")) # complex
    } else {
      form.glm.c=as.formula(paste(paste("Presence~",paste(paste0("poly(",names(env.stk_i),",3)"),collapse="+")))) # complex
    }
    
    # GAMs
    form.gam.s=as.formula(paste("Presence~",paste(paste0("s(",names(env.stk_i),",df=2)"),collapse="+"))) # simple
    form.gam.i=as.formula(paste("Presence~",paste(paste0("s(",names(env.stk_i),",df=3)"),collapse="+"))) # intermediate
    form.gam.c=as.formula(paste("Presence~",paste(paste0("s(",names(env.stk_i),",df=5)"),collapse="+"))) # complex
    
    # Trees
    form.trees=Presence ~ .
    
    ### define model settings
    modinp=list(multi("glm",list(formula=form.glm.s,family="binomial"),"glm-simple",step=FALSE,weight=FALSE),
                multi("glm",list(formula=form.glm.i,family="binomial"),"glm-interm",step=FALSE,weight=FALSE),
                multi("glm",list(formula=form.glm.c,family="binomial"),"glm-complex",step=FALSE,weight=FALSE),
                multi("gam",list(formula=form.gam.s,family="binomial"),"gam-simple",step=FALSE,weight=FALSE),
                multi("gam",list(formula=form.gam.i,family="binomial"),"gam-interm",step=FALSE,weight=FALSE),
                multi("gam",list(formula=form.gam.c,family="binomial"),"gam-complex",step=FALSE,weight=FALSE),
                multi("gbm",list(formula=form.trees,
                                 distribution = "bernoulli",
                                 interaction.depth = 5,
                                 shrinkage=.005,
                                 n.trees = 500),"gbm-simple",weight=FALSE),
                multi("gbm",list(formula=form.trees,
                                 distribution = "bernoulli",
                                 interaction.depth = 5,
                                 shrinkage=.005,
                                 n.trees = 1000),"gbm-interm",weight=FALSE),
                multi("gbm",list(formula=form.trees,
                                 distribution = "bernoulli",
                                 interaction.depth = 5,
                                 shrinkage=.005,
                                 n.trees = 10000),"gbm-complex",weight=FALSE),
                multi("randomForest",list(formula=form.trees,ntree=500,nodesize=10),"rf-simple"),
                multi("randomForest",list(formula=form.trees,ntree=500,nodesize=3),"rf-interm"),
                multi("randomForest",list(formula=form.trees,ntree=500,nodesize=1),"rf-complex"))
    
    ### =========================================================================
    ### Run models
    ### =========================================================================
    
    modis=wsl.flex(x=pseu.abs_i,
                   replicatetype="block-cv",
                   reps=3,
                   strata=blks,
                   project="tree_map",
                   mod_args=modinp)
    
    ### =========================================================================
    ### evaluate
    ### =========================================================================
    
    evals<-wsl.evaluate(modis,crit="maxTSS", prevalence_correction =T)
    
    smev=summary(evals)
    ord=sort(smev["tss",],decreasing=T)
    top6=which(colnames(smev)%in%names(ord)[1:6])
    modinp_top=modinp[top6]
    
    ### =========================================================================
    ### fit models for prediction
    ### =========================================================================
    
    prmod=wsl.flex(x=pseu.abs_i,
                   replicatetype="none",
                   reps=1,
                   project="tree_map",
                   mod_args=modinp_top)
    
    
    ### =========================================================================
    ### Predict
    ### =========================================================================
    
    ### Get thresholds
    thrs=get_thres(evals)
    
    # Null prediction to assess change
    prenull=wsl.predict(prmod,
                        predat=env.stk_i,
                        thres=thrs[top6],
                        write=TRUE,
                        output_dir=file.path(output_dir,"species_prediction"))
    
    ### =========================================================================
    ### combine
    ### =========================================================================
    
    # Load saved binary predictions
    if (!dir.exists(file.path(output_dir,"species_prediction"))) {dir.create(file.path(output_dir,"species_prediction"), recursive = T)}
    flbin=list.files(file.path(output_dir,"species_prediction",spi_name),full.names = T,pattern="Predictions")
    stk=stack(flbin)
    
    # Combine and save
    if (!dir.exists(file.path(output_dir,"ComitteeVote"))) {dir.create(file.path(output_dir,"ComitteeVote"), recursive = T)}
    flcom=file.path(output_dir,"ComitteeVote", spi_name, ".tif")
    calc(stk,sum,filename=flcom,datatype="INT1U",overwrite=TRUE)
    
    print(paste0("COMPLETED: [ ",round(id/length(Spname),2)*100, "% ]"))
    tmp <- tempfile()
    do.call(file.remove, list(list.files("tmp", full.names = TRUE)))
  }}

unlink(file.path(output_dir, tempfiles), recursive=T)
